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1.  Introduction 


Networks  of  finite  element  beam  elements  are  being  explored  as  a  means  of 
representing  trabecular  bone  structure  in  an  adaptive  hierarchical  multiscale  (HMS) 
simulation  framework.  The  beams  are  intended  to  comprise  a  representative 
volume  element  (RVE)  at  the  lower  length  scale,  which  is  interrogated  to  produce 
a  response  used  by  a  discretized  skeletal  model  at  the  higher  length  scale.  Use  of  a 
beam  network  will  facilitate  representation  of  an  anisotropic  structure  and 
variations  in  bone  density  by  modifying  the  beam  arrangement  and  cross-sectional 
properties.  The  hierarchical  multiscale  framework  is  described  in  detail  by  Knap 
et  al.1 

In  the  hierarchical  multiscale  framework,  the  models  functioning  at  both  the  higher 
and  lower  length  scale  exchange  information.  For  the  coupled  finite  element  models 
considered  here,  the  lower  length  scale  model  replaces  the  constitutive  evaluation 
in  the  higher  length  scale  code.  Consistent  with  this  treatment  as  a  constitutive 
model,  the  higher  length  scale  model  provides  a  set  of  history  variables  and  the 
velocity  gradient  to  drive  the  lower  length  scale  model,  and  the  lower  length  scale 
model  returns  Cauchy  stress  and  updated  values  of  the  history  variables.  A 
consistent  scheme  is  needed  for  specifying  boundary  conditions  on  the  RVE  and 
determination  of  the  average  Cauchy  stress  from  the  lower  length  scale  model. 
These  2  tasks  are  the  focus  of  this  report. 

In  the  construct  of  the  hierarchical  multiscale  framework,  the  lower  length  scale 
model  does  not  persist  from  one  time  step  to  the  next.  Their  persistence  would 
require  that  all  lower  length  scale  models  run  throughout  the  entire  simulation  and 
would  be  inconsistent  with  the  efficiency  objectives  of  the  hierarchical  multiscale 
framework.1  Instead,  for  each  time  step,  the  lower  length  scale  model  is 
reinstantiated  to  a  deformed  state  that  is  consistent  with  the  end  of  the  previous  time 
step.  The  reinstantiation  must  be  accomplished  using  a  reduced  set  of  information 
carried  by  the  higher  length  scale  model.  The  reduction  and  reinstantiation 
procedure  is  a  research  area  itself  for  the  many  applications  of  the  hierarchical 
multiscale  framework,  as  the  solution  is  often  path-dependent.  However,  for  a 
narrower  set  of  simulations  in  which  the  material  and  geometric  responses  are 
independent  of  deformation  path,  only  the  end  state  needs  to  be  specified.  The 
small-strain,  linear  elastic  examples  used  in  the  following  demonstration  problems 
are  assumed  to  be  path-independent. 

The  selection  of  an  appropriate  RVE  is  also  a  research  topic  that  will  not  be 
addressed  here.  The  term  “representative  volume  element”  is  used  in  the  most 
general  sense  of  representing  a  volume  of  material.  There  is  no  implication  that  the 
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homogenized  response  of  the  model  has  converged  with  respect  to  increasing  the 
volume  of  material  represented.  For  irregular  structures,  and  particularly  for  those 
with  gradients,  such  convergence  may  never  be  reached.  There  are  also  potential 
issues  that  arise  from  an  insufficient  number  of  beams  within  the  RVE  and  networks 
of  beams  whose  connectivity  results  in  otherwise  degenerate  structures.  Such 
instances  reflect  resolution  issues  and  model  construction,  and  those  types  of  sparse 
or  degenerate  beam  networks  are  not  considered  at  this  time. 

The  purpose  of  this  memorandum  is  to  document  and  validate  the  choices  made  for 
applying  boundary  conditions  to,  and  extracting  Cauchy  stress  from,  a  network  of 
elastic  beams  comprising  an  RVE  for  use  within  the  hierarchical  multiscale 
framework.1  While  it  has  been  shown  that  a  consistent  reduction  of  a  regular 
network  of  beams  can  lead  to  higher  order  gradient  models,2  the  determinations 
herein  are  restricted  to  simple  volume  averages  that  are  consistent  with  traditional 
finite  element  codes  that  may  be  used  as  the  higher  length  scale  model.  Higher  order 
representations  are  not  considered.  The  stress  averaging  outlined  in  the  following 
was  not  found  in  a  cursory  search  of  the  literature,  but  such  virtual  work  approaches 
are  standard  practice  and  there  is  no  assumption  that  it  is  original  or  unique. 

Vectors  and  second-rank  tensors  are  represented  by  bold  face  Latin  and  Greek 
characters,  respectively.  A  raised  dot  operator  indicates  an  inner  product,  and  a 
colon  operator  signifies  a  trace  of  the  inner  product  of  2  second-rank  tensors,  a 
contraction.  All  summations  are  explicitly  indicated,  and  no  summation  convention 
is  assumed. 

Unit  vectors  along  the  beam  directions  are  expressed  in  terms  of  the  coordinate 
direction  base  vectors  as  follows: 

—  "b  +  Tzez.  (1) 

S  Sx^X  T  SyBy  +  SZ6Z.  (2) 

t  tXBX  T  tyBy  +  tzBz.  (3) 

The  direction  of  the  beam  axis  is  denoted  by  the  unit  vector  r,  and  s  and  t  denote 

an  orthogonal  pair  of  vectors  that  are  also  orthogonal  to  the  beam  axis.  ex,  sy,  and 
bz  are  the  unit  vectors  along  the  global  Cartesian  coordinate  axes,  and  scalar 
projections  such  as  rx  —  r  ■  ex  are  direction  cosines  of  the  beam  unit  vectors. 

2.  Lower  Length  Scale  Model 

A  multitude  of  finite  element  codes  and  beam  elements  can  be  used  for  the  beam 
model,  and  the  boundary  conditions  and  stress  averaging  methods  described  here 
should  be  straightforward  to  apply  to  many  of  them.  The  Lawrence  Livermore 
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National  Laboratory’s  ParaDyn  code  is  used  within  the  HMS  framework,  so  it  is  an 
obvious  choice  for  the  present  evaluations.  ParaDyn  is  available  on  US  Department 
of  Defense  (DOD)  high-performance  computing  machines  without  restrictions  on 
the  number  of  simultaneous  executions.  It  has  fast  initialization  and  execution, 
robust  algorithms,  and  xdmf/hdf5  data  files  for  ease  of  extracting  requisite 
information. 

ParaDyn  is  an  explicit  dynamic  code  where  the  dynamic  stress  equilibrium 
equations  are  solved  by  a  second  order  accurate  time  integration  algorithm,  and 
information  is  propagated  through  the  structure  by  stress  waves.  The  time  step  is 
limited  by  the  Courant-Friedrichs-Lewy  condition  to  approximate  the  time  it  takes 
a  wave  to  traverse  an  element.  Hence,  many  small  time  steps  are  required  for  the 
solution  to  approach  a  quasi-static  equilibrium  configuration  as  waves  reverberate 
within  the  model  domain.  The  waves  and  dynamic  modes  persist  unless  quieted  by 
some  form  of  damping.  The  initialization  method  and  damping  described  in  the 
following  will  minimize  these  dynamic  effects  and  reduce  the  number  of  time  steps 
required  to  obtain  a  quiescent  solution. 

Beam  networks  are  generated  from  nodes  that  exist  either  on  the  exterior  boundary 
of  the  RVE  or  in  the  interior  of  the  RVE.  These  exterior  and  interior  nodes,  and 
their  connectivity,  form  the  beam  networks  of  the  lower  length  scale  model.  Unlike 
continuum  RVE  models  where  the  entire  exterior  of  the  RVE  is  populated  by  nodes, 
beam  structures  do  not  necessarily  define  a  nice  box  with  nodes  on  edges  and 
comers.  Beams  will  intersect  the  RVE  surface  at  multiple  angles  and  could  slide 
considerably  along  the  surface  if  constrained  only  in  the  direction  of  the  surface 
normal.  This  would  create  a  locally  high  average  deformation  on  the  surface 
inconsistent  with  the  applied  boundary  conditions.  Hence,  all  degrees  of  freedom 
of  all  exterior  nodes  are  constrained. 

The  specific  beam  element  used  in  the  evaluations  is  the  Hughes -Liu  beam,  which 
is  derived  from  a  consistent  reduction  of  a  shell  element.3  The  beam  supports  axial 
forces,  shear  forces,  bending  moments,  and  torsion.  The  beams  can  have  an 
arbitrary  cross-section  shape  (integrated  with  Gaussian  quadrature)  and  can  also 
have  a  linear  taper  along  their  length.  Only  circular  beams  of  constant  cross  section 
are  used  in  this  study.  The  method  presented  will  be  applicable  to  arbitrary  beam 
cross  sections,  but  additional  work  will  be  needed  to  determine  the  accuracy  of  the 
method  for  tapered  beams. 

Certain  assumptions  about  the  stress  distribution  are  built  into  the  beam  elements. 
While  stress  components  given  with  respect  to  the  beam  axes  represent  the  stress 
in  the  beam,  it  is  not  a  detailed  stress  field.  For  example,  the  parabolic  shear  stress 
distribution  across  the  beam  section  is  represented  by  an  average,  and  there  is  no 
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counterpart  to  this  shear  stress  resolved  along  the  beam  axis.  Consequently, 
standard  continuum  tensor  transformations  and  stress  averaging  do  not  apply 
because  of  assumptions  associated  with  the  dimensional  reduction.  This  is 
discussed  further  in  Section  4. 

3.  Reinstatiating  the  Beam  Network 

Within  the  HMS  framework,  a  reduced  set  of  information  is  supplied  to  the  lower 
length  scale  model.  This  information  is  used  to  first  reinstantiate  the  model  to  its 
prior  state  and  then  deform  it  to  the  current  state.  However,  this  process  can  result 
in  unstable  numerical  behaviors  if  the  path-dependence  is  important  and  not  taken 
into  account.  For  the  problems  considered  here,  linear  elastic  networks  of  beams 
are  subjected  to  small  deformations.  Elastic  solutions  are  path-independent 
provided  the  path  taken  does  not  encounter  a  material  or  geometrically  driven 
bifurcation  point,  such  as  buckling.  With  the  small  deformations  and  the  relatively 
small  beam  networks  considered  in  the  following  simulations,  no  bifurcation  points 
are  expected  in  the  vicinity  of  the  deformed  configuration,  so  the  final  state  is 
assumed  to  be  independent  of  deformation  path.  This  assumption  simplifies  the 
reinstatiation  problem  to  working  only  with  the  initial  beam  configuration  and  the 
current  deformation  specified  by  the  higher  length  scale  model. 

To  ensure  that  the  deformation  applied  to  the  model  is  consistent  with  the 
deformation  of  the  higher  length  scale  model  in  an  HMS  scheme,  all  of  the  nodes 
on  the  exterior  of  the  beam  network  model  are  moved  to  points  determined  by  their 
initial  coordinates  and  the  applied  deformation  gradient.  However,  abruptly  moving 
the  exterior  nodes  to  the  proper  locations  would  create  dynamic  effects  on  the 
interior  of  the  model  when  using  the  explicit  time  integration  scheme.  Alternatively, 
moving  the  exterior  nodes  slowly  enough  to  avoid  these  inertial  effects  would  be 
unacceptably  time  consuming.  The  approach  presented  here  is  to  move  both  the 
exterior  and  the  interior  nodes  to  positions  consistent  with  the  applied  deformation 
and  then  to  remove  the  constraint  from  these  interior  nodes.  This  allows  the 
coordinates  to  be  initialized  to  the  average  deformed  configuration  using  just  a  few 
time  steps,  and  the  subsequent  dynamic  solution  is  essentially  an  energy 
minimization  from  that  uniform  strain  state  with  the  exterior  nodes  fixed. 

Since  the  beams  are  elastic,  the  dynamic  motion  of  the  interior  will  persist  unless 
some  form  of  dissipation  is  applied.  Here,  a  viscous  force  is  applied  to  the  nodes 
through  a  mass  proportional  damping  algorithm  available  in  ParaDyn.  This  is  a 
nonconservative  force,  but  because  the  desired  quasi-static  solution  is 
path-independent,  the  precise  deformation  history  and  energy  loss  do  not  affect  the 
solution.  There  are  2  parameters  required  for  the  mass  proportional  damping  option: 
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the  peak  damping  frequency  and  the  fraction  of  critical  damping.  The  frequency 
can  be  estimated  by  observing  an  undamped  simulation  and  the  fraction  of  critical 
damping  is  chosen  to  be  0.9  to  allow  modest  overshoot  while  settling  into  the  static 
equilibrium  configuration. 

4.  Cauchy  Stress  Determination  from  Virtual  Work 

Several  approaches  can  be  used  to  determine  the  average  stress  over  the  RVE  beam 
model.  The  most  obvious  may  be  summing  forces  over  boundaries,  but  there  are 
potential  issues  because  the  deformation  is  applied  to  all  of  the  exterior  nodes.  The 
exterior  forces  in  the  3  coordinate  directions  will  sum  to  zero,  but  normal  forces  on 
opposite  faces  of  the  RVE  may  sum  to  different  values  because  of  tangential  forces 
on  the  remaining  surfaces.  This  does  not  produce  a  well-defined  average  stress. 
Similarly,  tangential  forces  on  corresponding  faces  may  not  balance  to  produce  a 
symmetric  shear  stress. 

Here  a  volume  averaging  method  based  on  virtual  work  will  be  employed.  It  is 
motivated  by  volume  averaging  stress  components  in  continuum  models,4  but  the 
implementation  details  are  somewhat  different  when  applied  over  beam  structures. 
The  method  assumes  that  an  equilibrium  solution  has  been  established  for  the  beam 
network  in  the  RVE.  The  average  stress  of  the  RVE,  a,  is  then  constructed  to  have 
the  same  virtual  work,  8<p ,  as  the  beam  network  when  subjected  to  a  common, 
work-conjugate  virtual  strain,  Se.  The  virtual  work  is  the  integral  of  mechanical 
work  over  the  RVE.  More  specifically,  it  is  the  volume  integral  of  the  contraction 
of  the  Cauchy  stress  tensor  with  the  work  conjugate  virtual  strain  tensor.  For  a 
constant  average  stress  over  the  RVE,  the  volume  integral  can  be  reduced  to  a 
simple  multiplication  by  the  RVE  volume. 

8<p  —  Jy  a\  8e  dv  =  a\  8e  V.  (4) 

Since  each  beam  is  a  discrete  element,  the  virtual  work  on  the  network  is 
determined  by  summing  the  contributions  from  the  individual  virtual  work  from 
each  beam  element. 


8(f)  =  Y"~Beams  f  0-Beam.  g£  dp  (5) 

vBeam 

Consistent  with  the  development  of  beam  models,  the  virtual  work  is  projected  into 
an  orthogonal  coordinate  system  aligned  with  the  beam,  and  the  volume  integration 
is  decomposed  into  integrals  over  the  beam  cross-sectional  area,  a,  and  the  beam 
length,  /. 
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50  =  El  Beams  fL  f  (ffA8eA  +  T s8ys  +  rT8yT )  da  d/,  (6) 

where  er^  is  the  stress  distribution  along  the  beam  axis  and  ts  and  tt  are  shear  stress 
distributions  over  the  beam  cross  section  along  the  local  orthogonal  directions  S 
and  T,  respectively.  The  virtual  strain  has  likewise  been  projected  in  this  local  beam 
coordinate  system  with  SsA  being  the  projection  along  the  beam  axis  and  8ys  and 
8yT  being  shear  strains  in  the  S  and  T  directions.  The  shear  strains  correspond  to 
gradients  along  the  beam  axis  of  virtual  displacements  in  the  S  and  T  directions. 
Shear  displacements  along  the  beam  axis  are  not  captured  by  beam  elements. 

The  beam  cross-sectional  area  is  assumed  to  be  constant  over  the  length  of  the 
beam,  which  allows  the  beam  length  to  be  integrated  separately.  The  integration  of 
the  stress  over  the  cross  section  is  accomplished  by  numerical  quadrature  for  the 
Hughes-Liu  beam  element  in  ParaDyn.  The  cross-section  integration  produces 
forces  in  the  beam  coordinate  direction  such  that  the  virtual  work  is  reduced  to  a 
summation  over  all  of  the  beams  in  the  RVE. 

50  =  If ~Beams  l  ( Fa8s  -I-  Fs8ys  +  FT8yT ).  (7) 

The  integrated  forces  FA,  Fs,  and  FT  appearing  in  Eq.  7  are  available  in  the  ParaDyn 
output. 

The  virtual  work  of  the  unknown  average  stress,  a,  represented  by  Eq.  4  is  set  equal 
to  the  virtual  work  given  by  the  explicit  summation  over  the  beams  in  Eq.  7.  Then, 
by  virtue  of  the  virtual  strains  being  arbitrary,  the  average  stress  can  be  determined. 
Determining  the  components  of  the  average  stress  requires  that  the  projections  of 
the  virtual  strain  tensor  appearing  in  Eq.  7  be  expressed  as  components  along  the 
global  Cartesian  coordinate  axes. 

The  virtual  strain  along  the  beam  axis  can  be  obtained  by  a  straightforward 
projection  of  the  strain  tensor.  The  definitions  of  the  projected  shear  strains  must 
be  consistent  with  the  virtual  work  expression  for  beam  elements,  Eq.  7.  The 
displacements  accompanying  the  virtual  shear  strains  are  orthogonal  to  the  beam 
axes  with  no  displacements  along  the  beam  axes.  This  is  reminiscent  of  simple 
shear  by  Sy.  Thus,  when  the  virtual  strain  is  projected  to  the  beam  coordinates,  only 
the  component  of  the  strain  associated  with  displacement  orthogonal  to  the  beam 
axis  contributes  to  the  virtual  work.  These  are 

8ys  —  s  ■  8e  r  and  SyT  =  t  ■  8e  r.  (8) 

The  virtual  strains  projected  along  the  beam  axis,  r,  and  in  the  5  and  t  directions 
are  expressed  as  components  along  the  global  Cartesian  coordinate  axes  using  the 
definitions  of  the  beam  vectors  from  the  notation  section. 
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(9) 


SeA  —  r  ■  8e  ■  r  —  rxex  ■  8e  ■  exrx  +  ryey  ■  8e  ■  eyry  +  rzez  ■  8e  ■  ezrz 
+2  ryey  ■  8e  ■  ezrz  +  2  rzez  ■  8e  ■  exrx  +  2  rxex  •  8e  ■  eyry 
?3c  xx  "E  9”  "E  2  TyTz8 £yz  +  2  YZTX8£ZX  "I"  2  TxTyS£xy. 

8y$  s  •  8e  •  v  vxsx8£xx  "I-  YySy8£yy  +  tzsz8£zz  -I-  (vysz-\-syvz')8£yz 

"EOz^je  "I”  SzT'x)^ £-zx  "E  ( J~x^y  "I”  ^x^yj^^xy.  (19) 

5 Yt  ~  t  '  ^E  '  V  —  Tx tx 8 £Xx  "E  Yyty8 £yy  +  Tztz8£zz  +  (vytz-\-tyVz^8  £yz 

"E( Tz^x  "l”  tzYx)8 £Zx  ”1”  (j"x^y  "1"  tx^y^S ^xy.  (11) 

Substituting  Eqs.  9-1 1  into  Eq.  7,  the  virtual  work  expression  for  the  beam  network 
is 

N-B  earns 

8<p  =  I  [(Fa^x  "E  Fsrxsx  +  Ft^x^x)  1  $^xx 
i 

+  {FAry  +  Fsrysy  +  Ftryty )  l  S£yy 
"E  (FArz  +  Fsrzsz  +  F{Vztz)  l  S£zz 

"E  (2  FAvyvz  +  Fsvysz  +  Fsvzsy  +  Ff:vytz  +  F^Tzty^  l  S£yz  (12) 
"E  (  2 Fatztx  +  Fsrxsz  +  Fsvzsx  +  F^vxtz  +  l  S£zx 

”E  (  2 FAvxvy  +  Fsvxsy  +  Fsvysx  +  Ff-Vxty  +  FfVyt-x^  l  S£xy. 

Eq.  12  is  simplified  by  introducing: 

Fx  —  Fatx  +  Fs  sx  +  Ft  tx 

Fy  —  FAry  +  Fs  sy  +  Ft  ty  (13) 

Fz  —  FArz  +  Fs  sz  +  Ft  tz  , 

so  that  the  expression  for  virtual  work,  Eq.  7,  becomes 

N-Beams 

8(f)  l  \FX^X^^XX  E  FyVySEyy  H~ 

1 

+(Fyrz  +  Fzry)8£yz  +  (Fzrx  +  Fxrz)8szx  +  (j Fxry  +  Fyrx)8exy\.  (14) 

Recalling  the  right-hand  side  of  Eq.  4,  the  virtual  work  of  the  average  stress  over 
the  RVE  subjected  to  the  same  virtual  strain  is  expanded  explicitly  as 

Sep  =  (t-.8eV 

1 7v:v 8 ('xx  E  oyy8 £yy  +  ozz8ez/  +  2  oyz8 £yz  +  2  ozx8  £/x  H-  2  oxy8 £xy]  E.  (15) 
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Since  Eqs.  14  and  15  are  equal  and  the  virtual  strains  are  arbitrary,  any  one  of  the 
strain  components  can  be  nonzero  while  setting  the  rest  to  zero.  The  stress 
components  for  the  RVE  can  then  be  determined  by  simply  matching  terms 
between  Eqs.  14  and  15.  Thus  the  average  Cauchy  stress  components  determined 
from  the  virtual  work  approach  are 


N-Beams 


OXX  y  ^ 


Fx  ^X  l  I 


N-Beams 

Oyy  y  Fy  Ty  l  f 

1 

N-Beams 

®zz  y  Fzrz  l  , 

1 

N-Beams 

Oyz  y  0.5  ( FyVZ  +  FZ  Ty  l  , 

1 

N-Beams 

®zx  ~  y  ^ '  0.5  ( Fzrx  +  FxVz)  l  ; 

i 

N-Beams 

0Xy  ~  y  ^  '  0.5  ( FXVy  +  FyVX ^  l  ■ 


(16) 


These  expressions  only  involve  the  forces  within  the  beams,  their  deformed  lengths 
and  orientations,  and  the  volume  of  the  RVE.  No  constitutive  assumptions  were 
required  for  this  derivation  implying  that  the  average  stress  calculation  is  applicable 
for  any  constitutive  model.  The  reinstantiation  procedure  for  these  cases,  however, 
would  have  to  be  revisited. 


5.  Verification  Analyses 

This  section  assesses  the  algorithm  for  determining  the  average  stress  within  the 
RVE  based  on  the  virtual  work  principle  (Eq.  16).  Since  the  averaging  process  itself 
is  the  topic  of  this  research,  the  best  available  comparisons  to  verify  that  it  produces 
favorable  results  are  comparison  with  either  the  total  elastic  energy  within  the  RVE 
or  the  summation  of  surfaces  forces  over  the  RVE  boundary.  Analytical  expressions 
are  derived  for  simple  cases;  however,  approximations  are  made  for  the  more 
complicated  examples.  It  is  important  to  note  that  inaccuracies  in  comparisons 
between  the  average  stress  calculations  and  the  validation  approximations  do  not 
necessarily  reflect  poorly  on  the  average  stress  calculation. 
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Rather,  as  will  be  demonstrated,  the  inaccuracies  reveal  limitations  of  the 
approximations  used  for  the  more  complex  geometries  and  therefore  a  limitation  in 
the  use  of  the  approximations  as  a  validation  tool. 

5.1  Beam  Properties  and  Loading  History 

The  beam  properties  and  sizes  used  in  the  analysis  were  chosen  for  convenience 
and  are  not  necessarily  representative  of  trabecular  bone  properties.  The  diameter 
of  the  round  cross  section  beams  is  0.1  mm,  and  the  Young’s  modulus  is 
E  =  400  MPa  (MPa  =  N/mm2).  The  yield  stress  is  set  to  100  MPa  so  that  the 
structure  remains  elastic.  Since  the  beam  deformation  is  integrated  through  time, 
the  logarithmic  strain  is  consistent  with  the  beam  deformation. 

In  all  cases,  the  peak  velocity  is  applied  to  the  RVE  for  1  ms  and  is  ramped  to  zero 
over  the  next  1  ms.  This  gives  a  total  displacement  in  millimeters  of  1.5  times  the 
velocity  in  millimeters  per  microsecond.  At  2  ms,  the  boundary  conditions  on  the 
interior  nodes  are  removed  and  the  interior  nodes  are  permitted  to  relax  to  their 
static  equilibrium  positions  over  the  next  2  ms.  In  the  following  examples,  the  peak 
velocities  are  all  less  than  10"13  mm/ms  at  the  end  of  the  4-ms  total  run  time.  This 
indicates  that  the  mass  proportional  damping  is  functioning  as  intended  to  eliminate 
the  dynamic  modes.  The  nodes  on  the  exterior  remain  at  the  locations  prescribed 
earlier  in  the  simulation. 

5.2  Analytic  and  Linearized  Energy  Relations 

Models  of  increasing  complexity  are  analyzed  to  assess  the  validity  of  the  average 
stress  calculation.  For  the  simplest  configurations,  the  stress  calculation  is 
compared  with  the  force  summation.  For  the  more  complex  beam  configurations, 
the  total  elastic  energy  of  the  beams  computed  by  ParaDyn  is  compared  with  the 
elastic  energy  integrated  over  the  RVE. 

e  =  f*[V(.e)<r(.e):de].  (17) 

The  volume,  V,  and  the  stress,  cr,  are  both  a  function  of  strain.  The  integration  limit, 
s,  is  the  average  RVE  strain. 

For  deformations  with  one  nonzero  strain  component,  only  a  single  conjugate  stress 
component  contributes  to  the  energy.  Thus,  if  the  applied  strain  and  energy  are 
known,  the  energy  can  be  used  to  validate  the  stress  calculation.  However,  the  stress 
is  not  necessarily  linear  in  strain  even  though  the  beams  themselves  are  linear 
elastic  and  the  strains  are  not  large.  The  geometry  changes  associated  with  the  beam 
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motion  introduce  nonlinearities.  Hence,  the  total  beam  energy  can  only  be  used  as 
an  approximate  verification  metric  for  the  virtual  stress  calculation  unless  the 
nonlinearties  of  the  system  can  be  properly  characterized. 

If  the  detailed  information  on  nonlinearities  is  unavailable  or  too  costly  to  process, 
a  linear  dependence  of  stress  on  applied  strain  is  typically  a  reasonable  assumption 
for  the  energy  in  a  linear  elastic  material.  For  uniaxial  strain  compression  where 
sev  is  the  applied  strain  and  oev  is  the  conjugate  stress  computed  by  Eq.  16,  the 
energy  is  calculated  by 


AL0e£d£j. 


(18) 


e  — 


The  exponential  times  the  reference  length  accounts  for  the  changing  volume  of  the 
RVE. 

5.3  Two-Beam  Model 

A  simple  2-dimensional  (2-D)  model  composed  of  2  beams  initially  oriented  at  45° 
with  respect  to  the  horizontal  is  constructed,  with  the  deformed  configuration 
shown  in  Fig.  1.  The  RVE  represented  by  the  bounding  box  is  H  (1  mm)  in  the  y- 
direction  and  W  =  2H  (2  mm)  in  the  v-di  recti  on.  Uniaxial  strain  velocity  boundary 
conditions  are  applied  to  force  compression  in  the  y-direction.  The  bottom  nodes 
are  fixed  and  the  2  upper  nodes  are  constrained  in  the  x-direction.  The  downward 
velocity  at  the  upper  surface  is  0.01  mm/ms  during  the  first  millisecond  so  that  the 
total  vertical  displacement  at  the  upper  nodes  after  the  ramp  down  in  velocity  is 
0.015  mm,  consistent  with  the  loading  description  in  Sections  3  and  5.1. 


Axiol  Force 


W=2H 


Fig.  1  Configuration  for  simple  2-beam  model  showing  the  axial  force.  H  is  the  initial  height 
(1  mm)  and  W  =  2H  is  the  initial  width  (2  mm). 
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5.3.1  Analytic  Solution 

The  deformed  beam  length  is 


l  = 


1 

(H  -  it)2  +  (y)  ]2  =  H  [(1  -  ft)2  +  1]2, 


(19) 


where  u  is  the  magnitude  of  the  downward  ^-displacement,  ft  is  the  displacement 
normalized  by  the  initial  height,  and  W  =  2H  has  been  utilized.  To  third  order  in  ft, 
the  log  strain  in  the  beams  is 

e  = -|(fi  “  (2°) 


The  force  in  the  vertical  direction  is  related  to  the  axial  stress  in  the  beams  ( aA  = 
E  s),  the  beam  area  (A),  and  the  current  angle  of  the  beam  with  respect  to  the 
horizontal  (9). 


Fy  =  oaA  sin(0)  —  E  e  A  — — - — T 

[(1-U)2+1]2 


-E  A  ^F(u--u2  -—u3  +  ■■■).  (21) 

2^2  V  2  24  ) 


With  an  initial  height  (H)  of  1  mm,  the  displacement  (u)  of  0.015  mm,  and  the  given 
modulus  and  beam  diameter,  Eq.  21  gives  a  force  of  0.016533824  N.  The  axial 
beam  force  computed  by  ParaDyn  is  shown  in  Fig.  1.  Multiplying  this  by  the  sine 
of  the  current  beam  angle  gives  the  same  force  as  Eq.  21  to  within  6  nonzero  digits. 
There  is  deviation  in  the  7th,  which  is  consistent  with  discarding  terms  0(u4)  in 
the  Eq.  21. 


The  average  y-direction  stress  on  the  RVE  is  the  force  in  y-di rcction,  multiplied  by 
2  for  the  2  beams  and  divided  by  2  mm  for  the  v-dimension  of  the  RVE.  A  unit 
depth  of  1  mm  is  assumed,  so  the  y-direction  stress  consistent  with  Eq.  21  is 
-0.016533824  MPa.  Analysis  of  the  beams  using  Eq.  16  gives  the  same  result,  with 
deviation  beginning  in  the  7th  nonzero  digit.  Similarly,  the  x-di rection  stress 
(determined  from  the  beam  force,  angle,  and  current  RVE  height)  is 
-0.017041248  MPa  and  the  shear  stress  is  0.016785629  MPa.  Both  agree  with  the 
results  of  Eq.  16  to  the  first  6  digits.  Hence,  the  stress  computation  using  the  virtual 
work  analysis  of  Eq.  16  is  accurate  for  this  simple  case. 


The  energy  per  unit  depth  for  the  2-D  model  is  obtained  for  each  beam  by 
integrating  the  force  times  the  differential  displacement.  The  result  is 

e  =  EA  -^=(-  u2  --u3  u4  +  (22) 

2V2  \2  6  96  )  ’ 


The  energy  computed  for  2  beams  from  the  analytic  energy  relation,  Eq.  22,  for  a 
displacement  of  0.015  mm  is  2.48647 4xl0'4  N-mm,  which  is  the  same  as  the  energy 
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integrated  by  ParaDyn  to  the  7  nonzero  digits  given.  This  establishes  that  a 
postanalysis  integration  of  the  energy  can  give  the  correct  result  if  all  of  the 
nonlinearities  are  considered. 

5.3.2  Approximate  Energy  Calculation 

However,  in  an  RVE  constructed  from  a  complex  beam  geometry,  there  will  not  be 
a  closed-form  solution,  and  it  is  important  to  know  the  magnitude  of  the  error 
introduced  by  not  directly  accounting  for  nonlinearities.  This  will  establish  how 
such  an  energy  metric  could  be  expected  to  be  in  error  when  verifying  whether  the 
virtual  work  stress  evaluated  from  Eq.  16  is  a  valid  means  of  computing  stress  for 
more  complex  geometries. 

For  the  previous  example  where  £ev  —  ln(l  —  0.015),  the  linear  energy 
approximation  of  Eq.  18  gives  e  —  0.0074811  A  L0aev  =  2.47383xl0"4  N-mm,  an 
error  of  0.5%.  A  simpler  approximation  of  e  —  0.5  A  L0aev  sev  =  2.49887x10  4 
N-mm  also  gives  an  error  of  0.5%. 

To  demonstrate  that  this  error  is  a  result  of  the  nonlinearities  in  the  problem,  the 
analysis  is  repeated  with  a  displacement  of  0.0015  mm,  an  order  of  magnitude 
smaller.  Here  the  y-di  reel  ion  stress  computed  by  the  virtual  work  method,  Eq.  16, 
is  -0.00166483  MPa  and  the  energy  integrated  numerically  by  ParaDyn  is 
2.49787xl0'6  N-mm.  Analytic  Eqs.  21  and  22  reproduce  these  results  to  at  least  6 
nonzero  digits.  However,  compared  with  the  results  of  the  0.015-mm  displacement, 
it  is  clear  that  the  stress  does  not  scale  linearly  with  displacement,  and  the  energy 
does  not  scale  quadratically.  The  energy  computed  by  the  linear  stress 
approximation  of  Eq.  18  is  2.49661xl0"6  N-mm,  an  error  of  0.025%,  which  is  an 
order  of  magnitude  smaller  than  with  the  larger  displacement.  These  results  are 
summarized  in  Table  1. 

Table  1  Average  RVE  stress  in  the  ^-direction  and  total  energy  computed  by  various  means 
for  the  configuration  shown  in  Fig.  1.  Digits  that  differ  between  stress  predictions  are 
highlighted.  Digits  of  the  energy  predictions  that  differ  from  the  ParaDyn  solution  are 
highlighted  in  yellow. 


Displacement 

Virtual  work 

Analytic 

Analytic 

Linearized  ParaDyn 

(mm) 

Eq.  16 

(t Tyy ,  MP3) 

Eq.  21 

((Tyy,  MPa) 

Eq.  22 
( e ,  N-mm) 

Eq.  20  ( e ,  N-mm) 

( e ,  N-mm) 

0.0015 

-0.00166483 

-0.00166483 

2.49787  |xl0'6 

2.49661xl0'6  2.497872xl0"6 

0.015 

-0.01653381 

-0.0165338| 

2.486474xl0"4 

2.47383x10  1  2.486474x10  ; 
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Hence,  at  the  0.015  strain  level,  errors  of  0.5%  in  calculating  the  total  energy  using 
the  linearized  approximation  of  Eq.  18  are  attributable  to  ignoring  nonlinearities.  If 
the  computed  energy  is  used  as  a  validation  metric  for  the  virtual  work  stress 
calculations,  agreement  to  within  0.5%  in  the  total  energy  is  the  best  that  can  be 
demanded. 


5.4  Irregular  Beam  Structure  in  Compression 


A  2-D  model  composed  of  111  beams  and  78  nodes,  Fig.  2,  is  used  for  further 
validation  of  the  virtual  work  stress  calculation.  The  model  extent  is 
2x2  mm  in  the  plane,  and  a  unit  length  out  of  plane  is  assumed.  Boundary 
conditions  consistent  with  uniaxial  compression  strain  are  applied  as  described  in 
Section  5.1,  with  the  total  displacement  on  the  upper  surface  being  0.015  mm.  The 
axial  stress  distribution  after  the  system  has  come  to  equilibrium  is  shown  in 
Fig.  2. 


Axiol  Force 
1.65e-03* 

-5.00e-03- 

-1.00e-02- 

-1.50e-02* 

-2-00e-02* 

-2  50e-02* 
-271e-02* 


Fig.  2  Irregular  beam  configuration  used  to  assess  the  stress  calculation.  The  axial  force  is 
depicted.  The  region  is  2  mm  in  both  the  x  and  y  directions. 

The  y-direction  stress  from  the  virtual  work  relation,  Eq.  16,  is  -0.04450714  MPa, 
and  the  full  stress  tensor  is  given  in  Table  2.  Substituting  this  stress  into  the 
linearized  energy  relation,  Eq.  18,  with  A  =  2  mm2,  Lo  =  2  mm  and  £ev  —  Zn(l  — 
0.015/2),  the  elastic  deformation  energy  is  estimated  as  6.66769x10  4  N-mm.  The 
integrated  elastic  energy  computed  by  ParaDyn  is  6.66622xl0‘4  N-mm,  an  energy 
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error  of  0.02%.  Therefore,  the  average  stress  for  a  complex  RVE  structure 
calculated  by  the  virtual  work  expressions,  Eq.  16,  is  consistent  with  the 
deformation  energy. 


Table  2  Computed  stress  for  irregular  beam  structure 


Stress  Components  for  j-direction  Strain  Stress  Components  for  Simple  Shear 

(MPa)  (MPa) 


-0.01675563 

0.00020244 

0.0 

0.00066804 

0.01038036 

0.0 

0.00020244 

-0.04450714 

0.0 

0.01038036 

-0.00004881 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Although  the  RVE  is  deformed  in  uniaxial  strain,  the  v-di  recti  on  stress  and  the 
xy-shear  stress  are  nonzero.  The  x-direction  stress  corresponds  to  the  stress  resisting 
Poisson  expansion.  The  shear  stress  results  from  the  reaction  of  the  structure  to  the 
applied  deformation.  The  structure  has  normal-shear  coupling,  a  natural  tendency 
to  shear  under  compression  loading.  However,  since  it  is  restrained  by  the  applied 
boundary  conditions,  a  reaction  stress  is  produced.  Neither  of  these  stresses 
contributes  to  the  macroscopic  energy  calculation  since  the  associated  strains  are 
zero. 


5.5  Irregular  Beam  Structure  in  Simple  Shear 

The  same  2-D  beam  model  shown  in  Fig.  2  is  subjected  to  simple  shear  boundary 
conditions  by  displacing  the  upper  surface  to  the  right  by  0.015  mm.  As  with  the 
compression  deformation,  all  of  the  nodes  are  first  moved  to  positions  consistent 
with  the  applied  strain  and  all  of  the  interior  nodes  are  subsequently  released  to 
relax  to  a  static  equilibrium  configuration. 

The  shear  stress  computed  from  the  virtual  work  relation,  Eq.  16,  is  0.01038036 
MPa.  Since  there  is  no  volume  change  associated  with  this  deformation,  the 
linearized  energy  relation  corresponding  to  Eq.  18  is  simply  e  —  0.5  A  L0aev  yev  = 
1.557054xl0~4  N-mm,  where  yev  —  0.015/2  is  the  shear  strain.  The  energy  directly 
integrated  by  ParaDyn  is  1 .557843x  104  N-mm,  for  an  error  of  0.05%.  Here  too, 
the  proposed  average  stress  calculation  is  shown  to  give  a  consistent  result. 

The  normal-shear  coupling  of  the  irregular  structure  is  also  evident  in  the 
x-direction  and  v-direction  stress  for  this  simulation,  shown  in  Table  2.  It  is  also 
notable  that  the  pressure,  -1/3  the  trace  of  the  stress,  is  not  zero,  implying  shear- 
pressure  coupling.  One  might  be  tempted  to  cross-check  the  2  analyses  of  the 
irregular  structure  by  calculating  the  K2212  entry  of  the  elastic  stiffness  from  the 
second  problem,  which  should  be  equal  to  the  K1222  entry  of  the  first  problem  if  the 
elastic  material  is  derivable  from  a  strain  energy  relation.  However,  simple  shear 
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involves  a  rotation  of  the  stress  tensor  to  satisfy  material  frame  indifference,  and 
that  contribution  to  the  x-  and  y-direction  stresses  is  expected  to  be  on  the  order  of 
approximately  10-4  MPa,  given  the  calculated  RVE  response.  This  unknown  stress 
contribution  is  similar  in  magnitude  to  the  stresses  involved  in  the  correlation,  so 
such  an  analysis  cannot  be  performed  with  any  degree  of  certainty  with  the 
information  currently  available. 

5.6  Deformation  in  the  Plastic  Range 

The  previous  analyses  demonstrated  that  the  average  stress  calculate  from  virtual 
work  relations,  Eq.  16,  provide  an  accurate  RVE  stress  in  the  elastic  range. 
However,  the  stress  calculation  does  not  depend  on  the  material  being  elastic,  as  it 
is  based  only  on  the  beam  forces  and  geometry.  The  stress  calculation  will  be 
equally  valid  if  the  beams  are  deformed  into  the  plastic  range.  This  is  clearly  evident 
for  the  2-beam  model  where  the  RVE  stress  is  computed  directly  from  the  forces 
and  geometry.  The  irregular  beam  structure  can  also  be  evaluated  in  the  plastic 
range,  but  the  energy  cannot  be  used  for  validation.  Instead,  validation  must  rely  on 
the  summation  of  nodal  forces  on  the  RVE  surfaces.  The  force  summation  requires 
that  the  nodes  on  the  lateral  surfaces  be  free  of  shear  tractions  (Section  3),  which  is 
only  appropriate  for  a  limited  range  of  deformation  modes,  as  will  be  discussed 
Section  6. 

A  simulation  similar  to  that  of  Section  5.4  was  run  with  shear- free  lateral  surfaces. 
If  the  beam  yield  strength  remains  high,  as  it  was  for  the  analyses  above,  the 
response  is  elastic  and  the  y-direction  stress  determined  from  the  virtual  stress 
relations,  Eq.  16,  is  0.03866  MPa.  The  force  output  in  ParaDyn  plots  is  available  to 
only  3  digits.  The  sum  of  the  y-direction  forces  on  the  upper  surface  divided  by  the 
2-mm2  surface  area  agrees  with  this  stress  to  the  first  3  digits,  as  expected.  The 
analysis  was  repeated  with  a  beam  yield  strength  of  0.5  MPa,  and  the  structure 
deformed  plastically.  Here  the  y-direction  stress  computed  from  the  virtual  work 
expressions  is  0.004542  MPa.  The  stress  computed  from  the  summation  of  surface 
forces  again  agrees  to  the  3  available  nonzero  digits,  confirming  that  the  method 
does  indeed  calculate  the  correct  stress  in  the  plastic  range. 

The  mass  proportional  damping  applied  to  quell  the  dynamic  modes  alters  the  path 
dependence  of  the  solution  in  the  plastic  range.  While  the  path  dependence  of  the 
solution  is  incorrect,  the  resulting  configuration  is  still  a  valid  static  equilibrium 
configuration,  and  it  is  suitable  for  analysis  stress  averaging  methods  for  the  RVE. 


15 


6.  Discussion 


The  decision  to  impose  the  boundary  conditions  on  all  of  the  exterior  nodes  of  the 
RVE  requires  discussion.  For  uniaxial  strain,  allowing  the  nodes  freedom  to  slide 
along  the  RVE  surfaces  does  result  in  somewhat  lower  stress  and  elastic 
deformation  energy,  as  noted  in  Section  5.6.  These  boundary  conditions  were 
investigated  and  the  results  are  also  consistent  with  stresses  calculated  by  the  virtual 
work  expression,  Eq.  16.  However,  it  is  not  possible  to  apply  shear  free  boundary 
conditions  consistent  with  an  arbitrary  strain  increment,  as  shear  stresses  cannot  be 
supported,  and  a  simple  shear  analysis  would  not  be  possible.  Thus,  an  arbitrary 
strain  can  only  be  imposed  consistently  on  the  RVE  by  specifying  boundary 
conditions  on  all  of  the  degrees  of  freedom  on  all  of  the  exterior  nodes. 

A  consequence  of  these  boundary  conditions  is  that  the  stress  cannot  be  calculated 
by  summing  forces  over  faces  except,  in  special  circumstances.  Consider  the 
configuration  in  Fig.  2.  The  y-direction  forces  on  the  lateral  surfaces  also  contribute 
to  the  y-direction  stress.  The  sums  of  the  forces  on  the  upper  and  lower  surfaces  are 
not  necessarily  the  same,  nor  is  the  sum  of  forces  over  a  face  divided  by  the  area 
necessarily  equal  to  the  average  stress. 

Another  boundary  condition  choice  that  needs  to  be  address  concerns  rotations  and 
moments.  No  rotational  constraints  are  imposed  on  the  beam  rotation  degrees  of 
freedom,  and  the  resultant  moments  at  static  equilibrium  are  zero  on  all  of  the 
nodes,  interior  and  exterior.  Imposing  rotation  constraints  on  the  exterior  of  the 
RVE  results  in  moments  on  the  exterior  of  the  RVE.  These  contribute  to  an  effective 
moment  imposed  on  the  RVE.  In  the  derivation  of  the  stress  tensor  for  continuum 
mechanics,  it  is  assumed  that  there  is  not  net  moment  at  a  material  point,  and  that 
assumption  results  in  the  stress  tensor  being  symmetric.  An  RVE  is  treated  as  a 
material  point  from  the  perspective  of  the  higher  length  scale  model,  so  having  a 
net  moment  on  the  RVE  would  be  inconsistent  with  the  underlying  continuum 
mechanics  assumptions. 

Simulations  on  the  beam  configuration  in  Fig.  2  were  run  imposing  constraints  on 
the  rotational  degrees  of  freedom  on  the  exterior  nodes  consistent  with  the  applied 
deformation.  For  simple  shear,  these  boundary  conditions  were  no  rotation  on  the 
upper  and  lower  surface  nodes,  and  nodal  rotations  equal  to  half  of  the  shear  angle 
on  the  lateral  surfaces.  The  rotations  of  the  comer  nodes  were  the  average  of  the  2 
intersecting  faces.  When  the  volume- average  stress  was  determined  using  virtual 
work  relations,  Eq.  16,  the  contribution  from  the  first  term  of  the  <rxy  expression 
was  different  from  the  contribution  from  the  second  term  in  the  first  digit.  Without 
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application  of  the  rotation  boundary  conditions,  the  terms  were  the  same  to  at  least 
7  nonzero  digits.  The  imposed  rotation  constraints  on  the  RVE  inconsistent  with 
the  underlying  continuum  mechanics  assumptions  and  the  virtual  work  derivation. 
The  classical  derivation  of  the  stress  tensor  is  for  a  material  point  and  assumes  that 
summation  of  moments  at  the  point  is  zero.  That  assumption  leads  to  the  stress 
tensor  being  symmetric.  Moments  resulting  from  imposition  of  rotation  boundary 
conditions  violate  this  underlying  assumption,  which  is  reflected  in  unequal 
contributions  to  the  shear  stress  in  Eq.  16.  Since  the  higher  length  scale  code  in  the 
HMS  framework  assumes  a  symmetric  stress  tensor,  rotational  constraints  cannot 
be  imposed  for  the  current  HMS  framework. 

Stepping  back  from  the  current  application,  there  is  nothing  inherently  wrong  with 
a  subscale  microstructure  model  producing  a  nonsymmetric  stress,  traction 
gradients  on  the  surfaces,  and  coupling  between  deformation  modes.  These  features 
reflect  aspects  of  microstructure  in  finite  volumes  that  are  not  considered  in  the 
continuum  mechanics  derivations  for  behavior  at  a  material  point.  The  upper  length 
scale  analyses  typically  ignore  the  structure  within  RVEs  because  only  the  first 
moments  of  the  response,  the  averages,  can  be  used  in  the  formulation.  It  should  be 
possible  to  construct  continuum- like  frameworks  with  finite- sized  base  units  rather 
than  assuming  reduction  to  a  material  point  (couple-stress  theories  are  still  based 
on  representations  at  a  material  point).  These  would  have  the  potential  to 
incorporate  second  moments  of  microstructure  information  into  large-scale 
deformation  analyses. 

7.  Conclusions 


An  efficient  procedure  has  been  presented  to  obtain  the  static  equilibrium 
configuration  of  a  network  of  elastic  beams  within  an  explicit  dynamic  finite 
element  code.  Expressions  for  computation  of  the  average  stress  over  the  RVE  have 
been  presented  based  on  the  principal  of  virtual  work,  and  these  stresses  have  been 
shown  to  be  consistent  with  both  nodal  force  computations  and  stored  energy.  The 
stress  averaging  procedure  is  also  valid  for  RVE  deformations  into  the  plastic  range. 
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